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^ ' We propose a theoretical model for ice growth under a wind-driven supercooled water 

film. The thickness and surface velocity of the water layer are variable by changing the 
Ch ■ 

^1 air stream velocity. For a given water supply rate, linear stability analysis is carried 



^ ■ out to study the morphological instability of the ice-water interface. In this model, 

water and air boundary layers are simultaneously disturbed due to the change in ice 

shape, and the effect of the interaction between air and water flows on the growth 

condition of the ice- water interface disturbance is taken into account. It is shown 

, Ph | that as wind speed increases, the amplification rate of the disturbance is significantly 

^ ■ affected by variable stresses exerted on the water-air interface by the air flow as well 
> ■ 

. as restoring forces due to gravity and surface tension. We predict that an ice pattern 

O 

O ' of a centimeter scale in wavelength appears and the wavelength decreases as wind 

cn ' 

^ ■ speed increases, and that the ice pattern moves in the direction opposite to the water 

^ ' flow. The effect of the air stress disturbance on the heat transfer coefficient at the 

water-air interface is also investigated for various wind speeds. 
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I. INTRODUCTION 



Thin liquid films are ubiquitous entities in a variety of settings and display interesting 
dynamics depending on various forces acting on the liquid film and the surface geometries on 
which the fluid moves.- A number of studies for the stability of a gravity-driven viscous liquid 
flowing on an inclined flat plane have been done, beginning from the pioneering works of 
Benjamin^ and Yih.- The stability of a wind-driven liquid film flowing over a horizontal flat 
plate^ and airfoils^i^ was investigated for small disturbances. In these works an interaction 
between air and liquid flows was considered, and it was shown that the thin liquid film 
becomes unstable to small disturbances, and that waves arise due to the variable stresses 
exerted on the liquid-air interface by the airflow. 

On the other hand, glaze (wet) ice formation and icicle growth are the problem of ice 
growth from a liquid film flow accompanying a phase change.- Glaze ice forms when water 
is collected from the impingement of supercooled water droplets, whereas icicles grow from 
the water of melting snow and ice at the root of the icicles. Typically, icicles also make up 
an important part of the total ice load in freezing rain.- The glaze ice and icicle surfaces 
are covered with a supercooled water film, and ice grows from a part of the water film by 
releasing latent heat into the ambient air below °C. It is well known that a solid surface 
under a supercooled liquid film is morphologically unstable, resulting in dendritic growth 
and a material is a microscopic mixture of solid and liquid. When a film of water is supplied 
by impinging droplets and cooled from its surface by cold air, the growing ice always initially 
entraps a considerable amount of liquid water. This is called spongy ice.-*^ It was recently 
shown that sponginess is a material parameter (70 % of ice), and is independent of the 
growth conditions.— The remaining unfrozen water flows on the ice surface. It should be 
noted that water flow is significantly different over an accreting ice layer than a non-accreting 
substrate and hence the problem of ice accretion in the presence of a flowing water film is 
highly complex. 

Ringlike ripples of a centimeter-scale in wavelength are formed.-^iiii^ Although the basic 
mechanism of icicle growth is well known,— >^ the mechanism of icicle ripple formation 
remained unsolved. Ogawa and Furukawa first attempted a theoretical explanation for the 
icicle ripple formation in the absence of airflow around icicles, where the icicle surface was 
covered with a gravity driven supercooled water film with a free surface.— In an improved 
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model proposed by Ueno,— — the influence of the shape of the water-air interface due to 
the action of gravity and surface tension on the ice growth conditions was newly taken into 
account, and a quite different mechanism on the origin of ripples on icicles was proposed. 
The theoretical results obtained by Ueno have recently been compared very favorably with 
experimental results.— We extended the above theoretical framework to include natural 
convection airflow around icicles.— It was found that the enhancement of the rate of latent 
heat loss from the water-air interface to the surrounding air due to airflow caused the 
amplification rate of the ice- water interface disturbance to increase. However, the wavelength 
of ice ripples was not significantly affected by the airflow. 

Since the natural convection airflow around icicles was less than 1 m/s in velocity, the 
free shear stress condition at the water-air interface was still satisfled, and hence driving 
force of the water fllm over the ice surface was gravity only.— However, when wind speed 
is large, the wind drag at the water-air interface also drives the water fllm. For example, 
the combination of two driving mechanism due to gravity and wind drag produces a variety 
of aufeis (also referred to as icings) morphologies with various surface features.— Aufeis are 
spreading and thickening ice accretions that form in cold air when a thin sheet of water 
flows or trickles over a cold surface. According to the aufeis formation experiments, an ini- 
tial morphologies of aufeis appeared essentially wavelike (or terraced) on a planar aluminum 
and a smooth ice surface, and their spacing and height varied with slope of an inclined 
plane and wind speed.— The morphological instability of the ice-water interface under the 
water fllm flow due to the two driving forces is also relevant to the surface roughness char- 
acteristics associated with glaze icing formation around aircraft wings and structures.— 
Furthermore, the morphological instability of the surface of growing ice is closely related to 
various natural phenomena where a thin layer of moving fluid separates the developing solid 
from the surrounding air.— 

The physical model commonly used in ice accretion codes is mainly based on conservation 
of energy and mass within numerical cells along the ice accreting surface.-i^ However, in 
glaze icing conditions the numerical results are poor agreement with experimental results. 
Therefore, some investigators developed theoretical and numerical understanding of the 
dynamic effect of water fllm flow on the glaze ice accretion. Bourgault et al. applied a simple 
fllm flow model to the problem of aircraft icing.— Myers et al. introduced a mathematical 
model for ice accretion with water flow driven by air shear, gravity, and surface tension, 
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employing the lubrication approximation to describe the water film flow.—"— A version of 
the Myers et al. model is used in the ICECREMO commercial aircraft icing code. The 
aerodynamic forces, as modified by the accreted ice, are significant in determining the wind 
drag and lift on iced structures. However, most analyses in current glaze ice accretion models 
lack the physical motivation for the effects of either surface roughness or profile change of 
ice on the heat transfer coefficient, and the roughness is treated as input to the code.— 
The lack of roughness formation in standard icing models indicates that the related surface 
instabilities must originate from more localized structures in which air flow can interact with 
the water film.— 

A more microscopic, rather than global, energy balance and detailed analysis of the in- 
teraction between the air and water flows are required to predict fine details of localized 
roughness.—"— Therefore, herein we perform a linear stability analysis for ice growth under 
a supercooled water film driven by a laminar airflow, taking into account the effect of inter- 
action between the air and water flows on the ice growth conditions. There is a significant 
difference between the current and previous works,—"— as follows: Our previous works fea- 
tured a gravity driven water flow, and the shape of the water-air interface was determined 
by the action of gravity and surface tension only. In the current model, water flow driven by 
air shear stress is considered. We will show that when the air and water flows are coupled, 
tangential and normal air shear stress disturbances as well as gravity and surface tension 
play an important role in determining the shape of the water-air interface as air stream 
velocity increases. Without employing any of the empirical methods used in standard icing 
models, the heat transfer coefficient at the water-air interface is determined explicitly by 
solving the governing equations for the air and water flows and the air temperature field. It 
will be shown that the growth conditions of the ice-water interface disturbance as well as 
the heat transfer coefficient at the water-air interface are strongly affected by the air shear 
stress disturbances, which is particularly a new effect not found in gravity driven water flow. 



II. MODEL 

The model configuration is shown in Fig. [H which is based on the experiments of aufeis 
formation on a smooth ice substrate in a wind tunnel by Streitz and Ettema.— In that 
experiment, water was supplied through a row of holes located at the upstream end of the 
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FIG. 1. Physical model of air- water-ice multi-phase system. Vertical height is not to scale. 

wind tunnel in a refrigerated laboratory, and the water film was driven by gravity and wind 
drag for various plane slopes and wind speeds. The following analysis is restricted to a 
two-dimensional vertical cross-section. The position x is measured from the leading edge 
where water is supplied at a constant rate, and the y axis is normal to it. The ice is covered 
with a thin water film oi Jiq. A cold air stream flows over the thin water layer, and the 
airflow is assumed to be laminar. The surface of the water layer moves at a velocity of uia 
under the influence of wind drag. The air velocity approaches the free stream velocity Uoo 
at a distance 6 from the water-air interface. Simultaneous water and air boundary layers 
occur. Since the air temperature at Too is lower than the ice-water interface temperature, 
Tsi, ice grows from a part of the water film by releasing the latent heat to the air through 
the water- air interface at temperature Tia- 

In this model, water is supplied from only the leading edge and there is no impingement 
of supercooled water droplets on the water film surface. The water and air boundary layers 
start at a; = 0. On the other hand, aircraft icing is primarily due to the impact of supercooled 
water droplets on a cold surface. One could see ice buildup at the leading edge of the airfoil, 
taken as a; = 0, where unfrozen water flow is slowest, but the water layer is thickest. In 
this sense, the model herein is not yet truly relevant to the aircraft icing problems.—"— In 
addition, the following assumptions are used in the current model: (1) The water film is 
driven by wind drag only. Hence the analysis is only valid on a horizontal surface, and 
the free stream velocity Mqo is constant in space. (2) Density remains constant through the 
phase change. (3) Change in ice shape disturbs the water-air interface, and the flow and 
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temperature fields in the water film and air. A quasi- stationary approximation is used for the 
disturbed fields and unsteadiness only enters through the Stephan condition, due to the long 
time scale of the ice-water interface motion. (4) Heat conduction into a substrate beneath 
ice sheet is not included. The ice sheet is assumed to be thick and the undisturbed part 
of temperature gradient in the ice does not exist. (5) The presence of waves on the water 
film is ignored because the waves did not interact with the forming ice in any observable 
manner in the experiments, except for enhancing the spreading of the water over the aufeis 
surface.— 



A. Governing equations 



The velocity components in the x and y directions in the air, Ua and w^, are governed 
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where pa is the air pressure, pa = 1.3 kg/m^, the density of air, and z/q = 1.3 x 10~^ m^/s, 
the kinematic viscosity of air. The velocity components in the x and y directions in the 
water layer, ui and vi, are governed by^ 

dui dui dui I dpi f d'^ui d'^ui\ 



dt dx dy pi dx \ dx^ dy 



dvi dvi dvi 1 dpi ( d'^vi d'^vi\ 

-^ + ui— + vi— = + ^i\^rT + ^TT ] - 9, (5) 



dt dx dy pi dy \ dx^ dy 

dui ^ dvi Q 
dx dy 

where vi = 1.8 x 10"^ m^/s and pi = 1.0 x 10"^ kg/m^ are the kinematic viscosity and density 
of water, respectively, pi is the water pressure and g the gravitational acceleration. The 
continuity equations ([3]) and ([6]) can be satisfied by introducing the stream functions ipa and 
il^i such that Ua = dtpa/dy, Va = -dtl^a/dx, ui = dtl^i/dy, and vi = -dtpi/dx. 



6 



Neglecting viscous dissipation in the energy equation, the equations for the temperatures 
in the air Tq, water T/ and ice Tg are^ 

dTg dTg dTg _ f d^Tg d^Tg \ 

dt dx dy " \ dx"^ dy"^ ) ' 

dTi dTi dTi _ ra^Ti d^Ti\ 
dt dx ^' dy '^^ \ dx"^ dy^ / ' 

where Kg = 1.87 x 10~^ m^/s, ki = 1.33 x 10"'' m^/s and Kg = 1.15 x 10"^ m^/s are the 
thermal diffusivities of air, water and ice, respectively. 



B. Boundary conditions 



The following boundary conditions are the same as those used in a previous paper- 
except for the first condition in Eq. ( !T9|) herein. Ignoring the density difference between ice 
and water, there is no normal fluid motion at the ice-water interface. Then both velocity 
components Ui and vi at a disturbed ice-water interface, y = ((t,x), must satisfy the no-slip 
condition:—"— 

Ul\y=(; = 0, Vl\y=^=0. (10) 

Since there is no impingement of supercooled water droplets on the water film, the kinematic 
condition at a disturbed water-air interface, y = ^{t,x), isr^'^ 



Vl\v=i- (11) 
The continuity of velocities of water film flow and airflow at the water-air interface is^ 

Ul\y=^ = Ug\y=^, Vl\y=i: = Vg\y=^. (l2) 

The continuity of tangential and normal stresses at the water-air interface leads to^'^'^ 
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where fii = pivi and Ha = Pa^^a are the viscosities of water and air, respectively, and 7 = 
7.6 X 10^^ N/m is the surface tension. The curvature term on the right hand side in Eq. 
(IT^ determines the magnitude of the surface tension induced stress. Hence, the condition 
expressed by Eq. (HM is that the capillary force resisting displacement and the normal stress 
on either side of the water-air interface should be in equilibrium.- 
The continuity condition of temperature at the ice-water interface is 

Ti\y=<: = Ts\y=(; = Ti, (15) 

in which the interfacial temperature Tj is an unknown to be determined. The Stephan 
condition is 



(16) 

s/=C 



dt J dy y=c dy 

which is based on the assumption that ice grows in proportion to the heat flux difference 
across the ice- water interface. Here L = 3.3 x 10^ J/m^ is the latent heat per unit volume, 
V is the undisturbed ice growth rate, and Kg = 2.22 J/(mKs) and Ki = 0.56 J/(mKs) are 
thermal conductivities of ice and water, respectively. 

The continuity condition of temperature at the water-air interface is 

^;|y=e = '^a\y=ii = Tla, (17) 

where Tia is the temperature at the water-air interface and will be determined later. The 
continuity of heat flux at the water-air interface is 

9r„ 



dy 



y=5 dy 



(18) 



where Ka = 0.024 J/(mKs) is the thermal conductivity of air. Far away from the air 
boundary layer, the velocities and temperature asymptote to their far-field values: 

'^a\y=oo ^oo; Va\y=oo 0, T(i\y=oo ^oo- (l'^) 



C. Linear stability analysis 

In this paper, the stability analysis will be limited to one-dimensional disturbances. A 
simple normal-mode analysis is applied to the field variables. Suppose an ice-water interface 
disturbance with a small amplitude (k resulting in ({ty^) = Cfc6xp[c7t + ikx], where k is the 
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wave number and cr = cr^'^^ + ia^^\ a^^^ and Vp = —a^^^'/k are the amplification rate and phase 
velocity of the disturbance, respectively. We separate ^, ipa, ipi, Pa, Pi, Ta, Ti and Tg into 
undisturbed steady fields with bar and disturbed parts with prime as follows: ^ = /iq + C 
i^a = i>a + ^'a, i^i = i^i + ^P'l, Pa = Pa + p'a. Pi = Pi + p'l, T„ = T„ + T^, = f , + T/ and 
Tg = Ts + T^. The undisturbed velocities in the air and water are derived from Ua = dipa/dy, 
Va = -dipa/dx, ui = dipi/dy, and vi = -dipi/dx. We define Ga = -dfa/dy\y=h^ and 
Gi = —dTi/dy\y=o as temperature gradients at the undisturbed water-air interface and ice- 
water interface, respectively. The disturbed field variables are assumed to be expanded in 
normal mode form, as follows: 



( 










Uoofaiv)^k 


i^'i 




Ulafl{y*)Ck 


p'a 




{paUlo/So)9a{v)^k 


p'l 




ipiul/ho)gi{y,)Ck 


K 




Ha{ri)Gaik 


V 




Hi{y.)GiCk 






\ Hs{y^)GiCk J 
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(20) 



We introduce the following two dimensionless variables tj = [y — ho)/6Q in the air and 
= y/ho in the water layer. Here Sq = {2i/ax/uooY^'^ is a scaled measure of the air boundary 
layer thickness,— u^o is the free stream velocity and x is the distance from the leading edge 
where water is supplied, uia in Eq. ( 120|) is the surface velocity of the water film driven by 
wind drag. As shown in Eq. f l29|l herein, Jiq and uia are functions of x. C,k is the amplitude 
of the water-air interface disturbance, and fa, fi, Qa, Qi, Ha, Hi and Hs are dimensionless 
amplitudes of disturbed parts of the stream function if), pressure p and temperature T. In 
the following, quasi- stationary approximation is used for the disturbed fields as in previous 
papers,—"— and we assume that the undisturbed part of temperature gradient within the 
ice does not exist, hence = T,/ (T^/ =0 °C for pure water). 
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1 . Equations and boundary conditions for undisturbed flows and 
temperatures in the air and water film 



Using the Blasius-type similarity transformations and substituting ipa = UooSoFa{ri) and 
Ta* = {Ta — Too)/ {Tia — T^) iuto the partial differential equations ([I]), ^ and ([7]), a set of 
ordinary differential equations for the dimensionless functions Fa and Ta* and their boundary 
conditions are obtained:— 
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(21) 
(22) 

0, (23) 



) Ta* I ri=0 1 ; 

where Pr^ = z/^/ is the Prandtl number of air. The first and second equations in Eq. f l23|) 
are derived from the undisturbed parts in Eq. (fT2|) by using the fact that the free stream 
velocity is much larger than the water surface velocity uia,-^ as shown in Table [B The 



third equation in Eq. ( 1231) is the result of the condition Ua\ 



y=oo 



dtpa/dy\y=oo = Uoo- The 



boundary conditions Ta 
Eq. m- 



y=ho 



Tia and Ta 



y=oo 



Too yield the fourth and fifth equations in 



For water film, we assume the following scaling = Cix", uia = C2X^ and T^; — T; 
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C^x'^, and ipi = uiah^Fi^y*). Here the constants Ci, C2, C3, a, h and c are determined from 
the boundary conditions, as follows. First, if the volumetric water flow rate per width. 
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(24) 



is constant, a + b = must hold, where ui* = Ui/uia = dFi/dy*. This is also derived by 
substituting ipi = uiahoFi{y*) into the undisturbed part of Eq. ( ITTi) . ui\y^j^^dhQ / dx = vi\y^i:^^. 
Second, the undisturbed part of Eq. ( |T3l) yields 
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■1/2 must hold. Finally, the undisturbed part of Eq. ( ITSl) yields 
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where Tia-T^o 
must hold. 



-Too is used because we assume \Tia\ <^ |Too|. From Eq. (126|) c — a = —1/2 
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Substituting ipi = uiahoFi{y*) and Ti* = (TJ — Tsi)/(Tsi — Tia) into the partial differential 
equations @ and ([8]), a set of differential equations for the dimensionless functions ui^, and 
Ti^, is obtained: 



Reiho _2 
2b— — —u 



d^Ti, Reiho _ ^ 
2ct^ — -FriUuTu, 



(27) 



(28) 



where i?e/ = Uiaho/ui and i^e^ = Moo^o/^^a are the Reynolds numbers of the water and air, 
and Pri = Ui/ni is the Prandtl number of water. Since ReJiQ/ [Rea^o) ^ 1 for values shown 
in Table m Eqs. (|27|) and ( p8|) can be approximated as d'^ui^/dyl = and d'^Ti^^/dyl = 0. 
The boundary conditions ui\y=o = 0, fiidui/dy\y=h^ = fiadua/dy\y=h^, T/|y=o = Tsi and 



T/q can be written as Mi*|y.=o = 0, dui^/dy.. 



*\y,=l 



1, Tj 



and T; 



/* I j/»=i 



— 1, respectively, by defining uia = HaUcxJi^d'^Fa/ drf\r^=Q/ {nido). Therefore, the solutions of 
the undisturbed velocity and temperature profiles in the water film are linear in y^, that is, 
Mj* = and Tl^, = — y*. This is in agreement with the more usual lubrication approach for 
describing icing with shear.-i^"— 

From a + b = 0, b — a = —1/2 and c — a = —1/2, we obtain a = 1/4, b = —1/4, 
c = —1/4. The value of a coincides with that stated in previous papers.-"^ Ci, C2 and C3 
are determined from Eqs. (124|) . (l25l) and (l26l) . Hence ho and uia can be expressed as follows: 
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It is found that ho and depend on Q/l^ and Mqo, and vary slowly with x. We assume that 
the scaling of ho and uia for a; holds except for the very vicinity of water source. Equations 
(!2T|) and (!23|) yield d^Fa/ drf\r,=Q = OAT. The shear rate for the undisturbed water layer is 
then given by 
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(30) 



y.=0 \2VaX ) III drf 

and its value is in the range 30.6 to 449.3 s~^ for the range of u^o = 5 to 30 m/s at x = 0.1 
m. Hence the value of the time defined by the inverse of shear rate lies within the range of 
2.2 X 10^'^ to 3.3 X 10~^ s for the above parameters. 
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Using Eq. fl2B]) . Ci = h^x and Sq = (2z/ax/uoo)"^'^^, we can express 

Tsi-T,^ = -C^^ G^^T^x-'/' = -^J^T^, (31) 

where Ga* = —dTa^/dr]\rj=o. Since the undisturbed part of temperature gradient within 
the ice does not exist in the model herein, the undisturbed part of Eq. ( !T6|) yields LV = 
Ki(Tsi — Tia)/ho, into which Eq. ( l3T|) is substituted to obtain the undisturbed ice growth 
rate 

K T 

V = " °° . (32) 

H^o/Ga*) 

Equations ( 13 ip and (152]) are the same form as those in a previous paper,— but the value of 
Ga* is different. From Eqs. (Ell), (I22D and ([23]), we obtain Ga* = 0.413 for Pr^ = 0.7. The 
variation of Tia and V against Uoo is shown in Table HIl 



2. Equations and boundary conditions for disturbed flows and temperatures 
in the air and water film 



When the assumed forms of ipa and are substituted into the complete equations ([T]) , 
1 and ([7]), the differential equations for the amplitudes fa and Ha are obtained: 
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(34) 



where fca* = kdo is the dimensionless wave number normalized by the length 5q. The dis- 
turbed part of Eq. ( 1T2]) and the boundary conditions u'a\y=oo = (^'^'al ^v\y=oo = and 
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(35) 



»y=o ■ ' dr] 

We note that, as shown in Table [T], since the water surface velocity uia is significantly lower 
than the free stream velocity Uoo, ^a|y=5 = and Va\y=^ = are good approximation, from 
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which the first and second equations in Eq. ( 135|) are obtained. Furthermore, the disturbed 
part of Eq. f|T7|) and the boundary condition T^\y^^ = give 



Ha\rj=0 — 1, 



H. 



a I r]=oo 



0. 



(36) 



On the other hand, when the assumed forms of ipi and Ti are substituted into Eqs. (jl]) 
([5]) and ([8]) and neglecting the terms with Reiho/ {RcaSo) ^ 1, the disturbed parts yield the 
equation for the amplitudes fi and Hi: 
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(37) 



(38) 



where A;/* = kh^ is the dimensionless wave number normalized by the length /zq, Ti^.{y^) = — 
and Pei = Uiaho/Ki is the Peclet number. We note that Eqs. (137|) and (138|) are in the same 
form as found in previous papers,—""— but the form of u/* is different. In this paper, ui^ = y^ 
is used. Using Eq. (129|) . the Reynolds number and the Peclet number of the water can be 
written as Rei = {2/ vi)Q /l^ and Pei = (2/ki)Q/Iw, which are independent of Mqo and x. 

Linearization of the disturbed parts of Eq. ( ITOl) at y = and Eqs. (fT3l) and ( |T4l) at 
y = ho yield the boundary conditions for ff. 



dfi_ 

dy* 

d'fi 



0, 



dyl 



y*=o 
d'fi 



//|y,=0 — 0, 



+ {kl + ^a)fl\y, = l=0, 



(39) 
(40) 



dyl 



- (Skl + ikuRei) 
y,=i dy^ 



+ikuRei 



Fr'' 



+ Wekt + Tla\fi\y^=i = Q 



where 



"d'fa 



Aii uia \So J dr]'^ 

Pa I ^oo 



Pi \UlaJ 5o 1 + {ka*ReaY \ df] 



ri=0 



r?=0 



(41) 



(42) 



3k 



2 " 



r]=0 



(43) 



Fr = uia/ {ghoY/"^ is the Froude number, and We = l/{piufJio) is the Weber number. 
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Linearization of the disturbed part of Eqs. ( !T7|) and (fT8|) at y = ho yields 



Hl\y^ = l + fl\y^ = , = 0, (44) 



dHi 



ho I dHa 



ri=0 



fiL=i = 0. (45) 



y,=i 6o \ dr] 

In deriving Eqs. ( 140|) . ( 14T]) . ( 144|) and ( l45l) . the relation between the amplitude of the water- 
air interface and that of the ice-water interface, = —fi\yt=iCk, is used, which is derived 
from the linearization of Eq. (fTT]) at y = /iq-— It should be noted that the water flow 
is affected by the terms in PUI) and Ua in (HTl) due to the tangential and normal air 
shear stress disturbances, respectively. The coupling between the air and water flows affects 
the disturbed part of temperature distribution in the water layer through the boundary 
conditions found in Eqs. (jS]) and (H5l) . 



3. Dispersion relation 

The disturbed parts of Eqs. (fT5|) and (fT6l) give the dimensionless amplification rate 
ai^^ = a^'^y {KaGal LJiq), and the dimensionless phase velocity fp* = —a^^'>/{kKaGa/L),—'— 



dy^ 



+ Kfki,{Hl'^\y,=o-l), (46) 




+ i^fA;,,i7«|,.=o), (47) 

y,=0 



where H^^^ and Hi^^ are the real and imaginary parts of Hi, and K[ = Kg/Ki = 3.96 is the 
ratio of the thermal conductivity of ice to that of water. 

The numerical procedure for calculating Eqs. (146|) and (147|) is as follows: First, Eqs. 
f l2T|) . fl22|) . f l33|) and are simultaneously solved for a given Uqo and x with boundary 
conditions ( 12^ . and The derived solutions Fa and are substituted into Eqs. 

dig) and (iSD. Using the boundary conditions ([3SD, (HHD and dH]), Eq. (E?]) is solved. Then 
Eq. ( 155]) is solved with the boundary conditions and (H^ . using solutions // and Ha- 
Finally, substituting solution Hi into Eqs. (146|) and ( 147|) and replacing ki^, with (ho/6o)ka*, 
Eqs. (j46|) and ( H7|l are obtained with respect to ka*. 
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TABLE I. Variation of a length, 5o, thickness of water film, Ilq, water-air surface velocity, uia, 
inverse of square of the Froude number, = gho/uf^, the Weber number, We = "y / (piufjio) , 

the Reynolds number of air, Rea = UooSo/i-'a, against the free stream velocity, Uoo, for Q/lw = 1000 
[(ml/h)/cm] and x = 0.1 m. The corresponding values of the Reynolds number and the Peclet 
number of water are Rei = uioIiq/vi = 31 and Pei = uiahQ/m = 418, respectively. 



(m/s) 


6o (/im) 


ho (/im) 


uia (cm/s) 


1/Fr2 


We 


Rea 


5 


721 


1348 


4.1 


7.78 


33.19 


277 


10 


510 


802 


6.9 


1.64 


19.74 


392 


15 


416 


591 


9.4 


0.66 


14.56 


480 


20 


361 


477 


11.7 


0.34 


11.74 


555 


25 


322 


403 


13.8 


0.21 


9.93 


620 


30 


294 


352 


15.8 


0.14 


8.66 


679 


100 


161 


143 


39.0 


0.01 


3.51 


1240 



III. RESULTS 

Figures [2] (a) and (b) show the variation of the values Wekf^, Tja \ Tia\ Ua^ and 

Ha'' against the water supply rate per width Q/lw in the range of 10 to 1000 [(ml/h)/cm], in 
the case of Uoo = 5 m/s and Uoo = 20 m/s, respectively. Here, T^a \ ni^''and T^a \ ni*^ are the 
real and imaginary parts of Eqs. f l42|) and f H3|) . respectively. From Eq. f l29|) . the variation 
of Eqs. (142|) . (143|) . and parameters Fr and We with respect to Q/l^ and follows: 
Sa ~ (Q/L)'/^ ~ (g//»)-^/^ l/Fr^ ~ iQ/Q-'^^ We ~ (Q/L)-'/^ for a given u^., 
and Tia ~ u^^^, Ha ~ u^^^, 1/Fr^ ~ u^^^, We ~ u'^^^ for a given Q/lu,- It is found 
that as Q/lw increases, the value of Sa increases, while other parameters decrease. Also, 
as Uoo increases, the value of decreases much slower than the other parameters. Hence, 
Fig. [2] (b) shows that is not negligible compared to the other parameters as Q/lw and 
increase. Therefore, we use the value Q/lw = 1000 [(ml/h)/cm] throughout this paper, 
which is also in the same order as that employed in the experiments.— 

On the other hand. Fig. [2](c) and (d) show the variation of the values 1/Fr^, Wekf^, 
T^a \ and Ha"* against the dimensionless wave number ka* for Q/lw = 1000 [(ml/h)/cm]. 
When plotting Wekf^ with respect to /ca*, the relation ku = {ho/8o)ka* is used. In the case 
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of Moo = 5 m/s, as shown in Fig. |2] (c), l/Fr"^ and Wekf^ are dominant terms in Eq. (HTI) . 
As Uoo increases, the values of l/Fr"^ and We decrease as shown in Table [H For example, for 
Moo = 20 m/s. Fig. [2] (d) shows that the magnitude of and Ua in Eqs. (HUj) and (HT]) are 
not negligible compared to l/Fr"^ and Wekf^. In Figs. [31 IllandElwe consider two cases: One 
takes into account the effect of the air stress disturbances and Ua, and the other does not. 
We will demonstrate in the following sections that this leads to critically different results 
for the shape of the water-air interface, the growth conditions of the ice-water interface and 
the heat transfer coefficient at the water-air interface. 




FIG. 2. Variation of l/Fr^, Wekl, T^P , ni''^ and ni*^ against Q/l^ for ka* = 0.2 and x = 0.1 
m, in the case of (a) Uoo = 5 m/s and (b) Uqo = 20 m/s. The variation of 1/Fr^, Wekf^, si''^ 
Ha ^ and Ha^ against /ca* for Q/lw = 1000 [(ml/h)/cm] and x = 0.1 m, in the case of (c) Uoo = 5 
m/s and (d) Uoo = 20 m/s. Here /cq* = 0.2 corresponds to a wavelength of about 2 cm for Uoo = 5 
m/s and about 1 cm for Uoo = 20 m/s at x = 0.1 m. 
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A. The shape of the water-air interface 




FIG. 3. For Q/lw = 1000 [(ml/h)/cm], Uoo = 5 m/s and Uoo = 20 m/s, (a) represents the variation 
of amplitude against ka*. Here /ca* = 1-0 corresponds to a wavelength of about 4.5 mm 

for lioo = 5 m/s, and about 2.3 mm for Uoo = 20 m/s at x = 0.1 m. (b) represents the variation of 
phase difference 0^^ between the ice- water and water-air interfaces against the free stream velocity 
Uoo for the most unstable mode. The solid curves consider the effect of the tangential and normal 
air shear stress disturbances on the water-air interface. The dashed curves do not consider this 
effect. 

It is supposed that a dimensionless small disturbance of the ice-water interface has a 
sinusoidal form: 

y* = C* = 5blm[exp{aJ^ + ik^x^)] = sm[ku{x^ - VpJ^)], (48) 

where 6b = Cfc/^o is an infinitesimal initial amplitude, a* = a/ (KaGa/Lho), = (V /ho)t, 
X* = x/Jiq, = 6bexp{ai^h^,) , and Im denotes the imaginary part of its argument. Since 

the water film is very thin, the deformed ice-water interface causes a disturbance of the 
water-air interface: 

= ^* = 1 + Im[5texp((T*t^, + iki^:X^)] 

= 1 - 6b{t^) sm[ku{x^ - VpJ^)] + //'^|y.=i cos[ku{x^ - v^*)]} 

= 1 + Sb{t^)\fi\y,=i\ sm[ku{x^ - Vp^U) - Q^X (49) 
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where/; and/; are the real and imaginary parts of/;, respectively, |/;|j^_^=i| = [(/; |?/,=i)^+ 
(//'-* |j;,=i)T''^ is the amplitude and cos e^, = -//''^|j^^=i/|//|j^^=i|, sinG^^ = //*\^=i/|/;|y,=i|, 
0^^ represents a phase difference between the water-air and ice-water interfaces. When de- 
riving the second equation in Eq. ( I49l) . the relation 5t = C,k/ho = —fi\y,=iSb is used.— 

Since fi\y,=i depends on the wave number, the amplitude and phase change according to 
the wavelength of the ice- water interface disturbance. Figure [3] (a) shows that the water-air 
interface tends to become flat as ka* increases, due to the action of gravity, surface tension 
and tangential and normal air shear stress disturbances on the water-air interface. In the case 
of Moo = 5 m/s, since the effect of air shear stress disturbances can be neglected, as shown 
in Fig. [2] (c), gravity and surface tension are dominant resisting forces for displacement of 
the water-air interface. On the other hand, as Uoo increases, a region of ka^, appears, where 
the action of tangential and normal air shear stress disturbances on the water-air interface 
is dominant compared to that of gravity and surface tension, as shown in Fig. [2] (d). For 
example, in the case of Mqo = 20 m/s in Fig. |3] (a), there is a region in the solid curve where 
the amplitude does not decrease with an increase in ka*- If we neglect the air shear stress 
disturbances, the amplitude is overestimated as shown by the dashed curve in Fig. [3] (a). 
As ka* increases, the difference between the solid and dashed curves becomes small because 
the surface tension Wekf^ is finally most dominant. 

Figure [3] (b) shows the variation of the phase difference B^^ between the ice-water and 
water-air interfaces against u^o for the most unstable mode (see IIII Bl) . with (solid curve) 
and without (dashed curve) the effect of air shear stress disturbances. In the case of Uoo = 5 
m/s, an upward phase shift of the water-air interface relative to the ice- water interface is 
large, as shown in Fig. |5](a). The solid curve shows that the phase difference decreases with 
increasing u^o and the sign of O^, changes from negative to positive at about u^o = 27 m/s. 
An example of the configuration of two interfaces at Uoo = 20 m/s appears in Fig. |5] (b). 
The decrease of phase difference B^, is also due to the effect of air shear stress disturbances 
on the water-air interface. On the other hand, if we neglect the effect of the air shear stress 
disturbances, the phase shift is still large even for large Mqo, as shown by the dashed curve 
in Fig. |3](b). 
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B. Amplification rate of the ice-water interface disturbance 



Figure 111(a) shows the variation of numerically obtained dimensionless amplification rates 
a^^'^ against the dimensionless wave number ka*- The solid curves represent (J^^\ taking into 
account the effect of tangential and normal air shear stress disturbances on the water-air 
interface. If we neglect the air shear stress disturbances, the dashed curves are obtained. 
In the case of Uqo = 5 m/s, the difference is negligible. However, if the air shear stress 
disturbances are neglected, the magnitude of a^J^ is overestimated with increasing Uoc- One 
expects to observe an ice pattern with a wave number at which the amplification rate is the 
maximum. For example, at Uoo = 20 m/s, a^"^ acquires a maximum value almax = 17.7 at 
ka* = 0.22. Since the wave number k is normalized by 6o, the corresponding wavelength 
of the ice pattern is 1.03 cm from A = 2Tx5Q/ka*. Here, the value of 5q = {2iyax/uooY^'^ = 
361 fim estimated from x = 0.1 m and Uoo = 20 m/s is used. The magnitude of Vp* is 
defined from the wave number at which ai^^ acquires a maximum value. At ka* = 0.22, 
we obtain Vp* = —96.6, and hence the displacement of the ice-water interface after the time 

= 1/crimax is Ax* = I'p^/crimax = —5.4. The icc pattern will move in the direction opposite 
to the water fiow (see Fig. [5]). The variation of crlmax, A, Vp* and Ax* against Uoo is shown 
in Table M 

It is found from Fig. H] (b) and Table [Til that the wavelength shortens with increasing 
Uoo- Wavelike ice patterns with various roughness spacings and heights were experimentally 
observed by changing the wind speed and slope of an inclined plane.— For a wind speed of 
16 km/h=4.4 m/s, the roughness spacing is increased as the plane slope is decreased (the 
roughness spacing for smooth-ice base is about 3 cm at about 3°, see Fig. 10 in Ref. 13) and 
for the plane slope of 8^, the roughness spacing is decreased as the wind speed is increased 



(see Fig. 11 in Ref. |20| ). The latter result is consistent with the theoretical prediction, 
except that the results herein are only those obtained with a slope of 0°. 

Since the values ai^"^ in Table IITI are much larger than those predicted by the morphological 
instability triggered by thermal diffusion at the water-air interface in previous papers,—"— 
the ice-water interface instability herein is enhanced by the flow in the water fllm. On the 
other hand, as shown in Figs. [2] (c) and (d), the surface tension Wekf^ is most dominant 
to suppress the water-air disturbance with increasing ka* and stabilizes the corresponding 
ice-water interface disturbance. The solid and dashed curves in Fig. H] (a) show that as 
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(a) 



/ 

/U„=10 m/s\ 




Uoo=5 m/s 



AOS 0.1 



Otis ai: 0,iS 0:3 




Uoo (m/s) 



FIG. 4. For Q/Z^ = 1000 [(ml/h)/cm] and x = 0.1 m, (a) dimensionless amplification rate cj* = 
a^"^^ / {KaGa/ LKq) versus dimensionless wave number ka* = k6o, (b) variation of wavelength of 
ripples against the free stream velocity Uoo- The solid curves consider the effect of the tangential 
and normal air shear stress disturbances on the water-air interface, and the dashed curves do not 
consider this effect. 

Uoo increases, the wave number at which ai^^ vanishes, that is, the neutral stability point is 
shifted to higher wave number. This is because the value of the Weber number We decreases 
with an increase in Uoo, as indicated in Table HI hence the stabilization due to the surface 
tension Wekf^ becomes more effective for higher wave numbers. As shown in Figs. |2] (c) 
and (d), since Ili^'' has negative values with respect to /Ca*, the value of Wekf^ + Ili'^'' in Eq. 
( HTl) decreases as u^o increases. This means that the stabilization due to the surface tension 
Wekf^ is weakened by normal air shear stress disturbance ni^'*. Hence, the wave number of 
the neutral stability point in the solid curves is shifted to the higher wave number compared 
to that in the dashed curves with an increase in Uoq. That is why the wavelength evaluated 
from the most unstable mode becomes shorter as u^o increases. It should be noted that the 

(r) 

magnitude of a* is decreased by the effect of the tangential and normal air shear stress 
disturbances. However, the effect of air shear stress disturbances on the wavelength does 
not make a significant difference, as shown in Fig. Il](b). 

On a static structure, typical wind speed is in the order of 20 m/s.— The wind speed in the 
aufeis formation experiments by Streitz and Ettema^^ was varied up to 48 km/h=13.3 m/s. 
On the other hand, in aircraft icing a typical value is around 100 m/s.— By applying our anal- 
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ysis to the case Mqo = 100 m/s, assuming that the air stream flow remains laminar as shown 
in Fig. [H the results shown in Tables |T] and [TTl are obtained. On the other hand, in the limit 
Uoo 0, there is no driving force to move the water film. In this case, Sq = {2vax/uooY^'^ 
and Jiq in Eq. f l2^ have an infinite value and hence the corresponding wavelength is not 
defined. The same issue arises for gravity driven water flows found in previous papers,—""— 
in which the thickness of water film is determined from Tiq = [Sui/^g sin 9)Q/lu]^^^ , where 9 
is the inclination angle with respect to the horizontal. In the limit 6* — )■ 0, there is no driving 
force to move the water film. Then the wavelength at 6' = is not defined (see Fig. 8 (a) in 



Ref. 
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TABLE II. Variation of temperature at the water-air interface, Ti^, undisturbed ice growth rate, 

- • (r) 

V, maximum value of dimensionless amplification rate (T*max at a dimensionless wave number ka*, 
the corresponding wavelength. A, dimensionless phase velocity, Wp*, dimensionless displacement of 
the ice- water interface, Ax^ = Vp*/crimax after the dimensionless time = l/uimax, against the 



free stream velocity, Uoo, 


for X = 0.1 m 


Q /Iw 


= 1000 [( 


ml/h) /cm] 


and Too = 


-10 °C. 




Uoo (m/s) 


Tia rc) 


V (xlO-^ 


m/s) 


k 


" *max 


A (cm) 




Ax* 


5 


-0.32 


0.40 




0.12 


18.1 


3.78 


-70.3 


-3.9 


10 


-0.27 


0.57 




0.17 


23.9 


1.88 


-97.9 


-4.1 


15 


-0.25 


0.70 




0.21 


21.0 


1.25 


-105.1 


-5.0 


20 


-0.23 


0.81 




0.22 


17.7 


1.03 


-96.6 


-5.4 


25 


-0.22 


0.91 




0.23 


16.0 


0.88 


-91.8 


-5.7 


30 


-0.21 


1.00 




0.23 


15.2 


0.80 


-80.1 


-5.3 


100 


-0.15 


1.83 




0.29 


13.8 


0.35 


-68.7 


-5.0 



C. Heat transfer at disturbed ice- water and water-air interfaces 

Using the assumed forms of T/ and in Eq. f l20|) and considering their imaginary parts, 
the temperature in the water layer and ice can be expressed in dimensionless form as follows: 

J- si — -L la 

+H^'\y*) cos[ku{x^ - Vp^U)]^ , (50) 
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J- si — J- la 

+Hl'^\y^=o cos[A;/*(x* - fp*t*)]| , (51) 

where we used T), = Tsi and the solution Hs{y^) = {Hi\y^=Q — l)exp(/c/*t/,,).— 

A microscopic energy balance have to be considered to explain fine details of ice mor- 
phology. We define the disturbed parts of heat fiux from the ice-water interface to the 
water, from the ice to the ice-water interface, and from the water-air interface to the air, 
as ql ^ lm[-KidTl/dy\y=c], q', ^ lm[- K ^dT^ dy\y=c] and q', ^ lm[-KadT:^/dy\y=^], respec- 
tively. These can be expressed in dimensionless form as follows: 



ii* 




sin[/c;^(X;f ^p*^*)] 

y,=0 



Kid I dy^ 

cos[ku{x^ — t'p*t*)] , (52) 



I's* = TTTT = -Sb{t*)Kiku \{H^''^\y^=o - 1) sm[ku{x^ - VpJ^)] 
J^i^i ^ 

+Hl'''^\y^=o cos[ki^{x^ - VpM)]^ , (53) 



q' = 



^^^^ = -SbiQ { (Cf V/^i.=i - //^^l..=i) sin[A;,,(x, - v^J^)] 

+ (Cf V/^i.=i + //'i.=i) cos[A;,.(x. - vpj,)]}, (54) 

where G'a'^ = {hQ/6o){—dHa^ /dri)\ri=o and Ga*-* = {ho/6o){—dHa^/dri)\rj=o represent the 
real and imaginary parts of the disturbed part of the air temperature gradient G'^ = 
{ho/6o){—dHa/dr])\ri=o at the water-air interface. From Eq. ( !T6ll . the disturbed part of the 
Stephan condition in dimensionless form can be written as = ql^ — q'^^. Substituting 

Eqs. (08]), (152]) and ([53]) into this condition, Eqs. (jM]) and (07]) are obtained. 

Figures |5] (a) and (b) illustrate the time evolution of the ice- water interface disturbance 
with an initial amplitude of 5b = 0.1 in the case of Uoo = 5 m/s and Uqo = 20 m/s, 
respectively. The respective wave numbers of disturbance are ka* = 0.12 in Fig. 151(a) and 
ka* = 0.22 in Fig. |5] (b). These are the fastest growing modes, at which ai^^ acquires a 
maximum value, as shown by the solid curves in Fig. HI (a). The phase shift of the water-air 
interface relative to the ice-water interface in Fig. [5] (b) is negligibly small compared to 
that in Fig. [5] (a), as shown by the solid curve in Fig. [3] (b). Due to the left-to right air 
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FIG. 5. For Q/lw = 1000 [(ml/h)/cm], (a) and (b) are illustrations of the time evolution of an 
initial disturbance of the ice- water interface from t^, = to = 1/ crlmax • The arrows indicate the 
position of maximum point of disturbed heat flux q'^^, q'^^ at the ice- water interface and that of q'^^ 
at the water-air interface, (a) represents the disturbance of ka^ = 0.12 in the case of Uoo = 5 m/s. 
(b) represents the disturbance of ka* = 0.22 in the case of Uoo = 20 m/s. Ax* is the displacement 
of the ice- water interface after the time t^, = 1 

^*max* Vertical height is not to scale. 

and water flows indicated by arrows in Fig. [5l the isotherms in the air and water boundary 
layers are no longer symmetrical around each protruded part of the water-air and ice-water 
interfaces. Since the isotherms become closer on the upstream side of each protruded part 
of the interfaces, g^^, and q'^^ are largest on the upstream side of each protruded part, 
as indicated by the vertical arrows in Fig. O Hence, the ice growth rate on the upstream 
side of each protruded part is faster than that on the downstream side, and this results in 
the translation of the ice-water interface in the direction opposite to the water flow. As 
mentioned in IIIIBt the displacements after the time = l/a^max are Ax* = —3.9 and 
Ax* = —5.4 in Figs. E] (a) and (b), respectively. 

We separate the local heat transfer coefficient at the water-air interface, hx, into 
the undisturbed part hx = —KadTa/dy\y^h^^/{Tia — Too) and the disturbed part h'^ = 
lm[-KadTjdy\y^-h^/{Tia-Too)\. The former can be written as hx = Ka/ 6o{-dfa*/dr]\jj=o) ■ 
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Using the value of Ga* = ~dTa*/dr]\r,=o = 0.413 obtained in III CI the value of local Nusselt 
number scaled by y/ReZ is Nu^/ VR^ax = {hxx/Ka)/ y/uoox/ua = -{1/ V2)dfa*/dr]\n=o = 
0.292, from which we obtain = 0.292Ka\/uoo/ (i^ax). A similar expression for the laminar 
convective heat transfer coefficient, hg = 0.296{Ka/ y/T^) ^ V^-^"^ / {J^ V^-^'^ds), is used in 
aircraft icing models,—*^ where s is the surface distance from the stagnation point and Ve is 
the velocity at edge of air boundary layer. When Ve = Mqo (constant) and replacing s with x, 
hg yields the same expression as except for the very slight difference of numerical factor. 
Using Tlx, the undisturbed ice growth rate V in Eq. (!32|) can be expressed as = —hxToo/L. 
On the other hand, the disturbed part of the heat transfer coefficient normalized by the 
undisturbed one can be written as 

K/hx = -Im [G'Ji\y^=i6bexp{a4^ + ik^x^)] 

= Sb{U) [{h'^/hx)'''''' sm[ku{x^ - v^J,^)\ + {KJhx)^'^ cosffc/^^^* - ^^p*4)]] 

= h{t^)\h!Jhx\ sinf/ci^x* - v^J.^) - O^/^J. (55) 

Equation (l55ll becomes the same form as Eq. ( 15^ by putting {h'^/lixY^''' = — (Ga'^V/'^'' |i/*=i ~ 
Gf^fl%=i) and {h'JhxY'^ = -{G'^''^f^%=^ + Gf^f^\,=^). Here 

\K/hx\ = [{{h'x/hxf^V + m/hxYnV (56) 

is the amplitude, and cos6g/^_^ = {h'^/hxY'^y\h'^/hx\ and sin6g/^_^ = —{h'^/hxY''''/\h'^/hx\, Qq'^^ 
is a phase difference between q'^^ = h'^/hx and the ice- water interface. 

The solid curves in Fig. [6] (a) show the variations of {h'^/hxY^^ and [h'^/hx]'^'^^ against 
Moo, taking into account the air shear stress disturbances. These values are estimated for 
Q/lw = 1000 [(ml/h)/cm] and x = 0.1 m, and for ka* at which ai^'' acquires a maximum 
value. It should be noted that q'^^ = h'^/hx includes G'^ and fi\y,=i. hx depends on only 
two parameters, free stream velocity u^o and position x. On the other hand, h'^ depends 
on many parameters. = (hQ/6o){—dHa/d'i])\ri=o is determined from the disturbed airflow 
and temperature flelds, and fi\y,=i determines the magnitude of amplitude and phase of the 
water-air interface by using the relation = —fi\yt=iCk- The shape of the water-air interface 
changes by the action of gravity, surface tension and air shear stress disturbances. As shown 
in Fig. [31 if we neglect the effect of the air shear stress disturbances, the amplitude and 
phase of the water-air interface relative to the ice- water interface are not correctly evaluated 
as Uoo increases. This results in an overestimated value of \h'^/hx\ compared to that taking 
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FIG. 6. For Q/lw = 1000 [(ml/h)/cm] and x = 0.1 m, (a) represents the variation of the disturbed 
part of the heat transfer coefficient normahzed by the undisturbed heat transfer coefficient against 
the free stream velocity Uo^'- (^z/^x)^^-* is real part and (h'^/hx)^^'^ is imaginary part, (b) represents 
the variation of \h'x/hx\ against Uqo- The sohd curves consider the effect of the tangential and 
normal air shear stress disturbances on the water-air interface, and the dashed curves do not 
consider this effect, (c) represents the variation of {h'^ /hx)^^'^ (solid curves) and {h'x /~hx)^'^^ (dashed 
curves) against dimensionless wave number ka* for free stream velocities ?ioo=5, 10, 20 m/s. 

into account the effect of air shear stress disturbances with increasing u^, as shown by the 
dashed and solid curves in Fig. |6] (b). 

The solid curves in Fig. [H](a) shows that (h'^/hx)^^^ is positive for any Uoo, while {h'x/hxY^'^ 
is negative when < 10 m/s and is positive when Mqo > 10 m/s. From cos6g/^^ = 
{h'x/hxY'^y\h'x/hx\ and sin6g/^^ = —{h'^/hx)^^y\h'x/hx\, the corresponding phase difference 
between q'^^ = h'x/hx and the ice-water interface is — vr < Qq'^^ < — 7r/2 and — 7r/2 < Qq'^^ < 0, 
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respectively. On the other hand, if we neglect the air shear stress disturbances, as shown 
by the dashed curves in Fig. |6] (a), {h'^/h^)^^^ is negative and (h'^/hx)^'^^ is positive for any 
Moo- Then, the phase difference between q'^^ = h'^/h^ and the ice-water interface is always 
— TT < ©g/^^ < —7i/2 with respect to any Uoq. This means that the position of maximum 
point of heat flux q'^^ or that of the heat transfer coefficient h'^/hx depend significantly on 
the air shear stress disturbances. 

Figure M (c) shows the variation of {h'^/hxY^'' (solid curves) and (h'x/hx)^^'' (dashed 
curves) against ka^, for the free stream velocities of Moo=5, 10, 20 m/s. In the case of 
Uoo = 5 m/s, {h'^/hxY^^ is negative and {h'^/hxY''^ is positive for any ka*. From cosQqi^^ = 
{h'^/hxY^'^ /\h'^/hx\ and sin6g^^ = —{h'^/hxY'''y\h'^/hx\, the phase difference between q'^^ = 
h'^/hx and the ice-water interface is — vr < 6^/^^ < — 7r/2 as shown in Fig. [5] (a). On the other 
hand, in the case of Moo= 10 and 20 m/s, {h'^/hxY^^ is positive for any ka*, but {h'^/hxY'^^ 
changes from positive to negative at ka^ = 0.17 and /Cq* = 0.33, respectively. Then, in 
the case of = 10 m/s, the phase difference is — 7r/2 < 0^/^^ < for ka* < 0.17 and 
—IT < ©g/^^ < — 7r/2 for ka* > 0.17. Likewise, in the case of u^o = 20 m/s, the phase differ- 
ence is — 7r/2 < Qqi^^ < for ka* < 0.33 and —n < 6^^^ < — 7r/2 for ka* > 0.33, as shown in 
Fig. [5] (b). This means that the position of maximum point of q'^^ = h'^/hx on the water-air 
interface changes according to the wavelength of the ice-water interface disturbance and the 
free stream velocity. Furthermore, as shown in Fig. \5[ the position of the maximum point 
of q'^^ = h'^/hx moves in the direction opposite to the water flow with time. 

IV. SUMMARY AND DISCUSSION 

We have proposed a theoretical model for ice growth under a supercooled water film 
driven by wind drag. The thickness and surface velocity of the water layer are variable 
by changing air stream velocity and water supply rate. For a given water supply rate, we 
investigated the morphological instability of the ice-water interface for various air stream 
velocities using a linear stability analysis, taking into account the effect of gravity, surface 
tension and the tangential and normal air shear stress disturbances due to the airflow on 
the shape of the water-air interface. 

Even for the simple model developed here, the form of heat transfer coefficient at the 
disturbed water-air interface is too complicated, which depends on the disturbed air flow 
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and temperature fields, the shape of the disturbed water-air interface, as well as the shape of 
the ice-water interface. By considering the interaction between the air and water flows, we 
have found that the heat transfer coefficient at the water-air interface is significantly affected 
by the air shear stress disturbances, which suppresses the dimensionless amplification rate 
of the ice-water interface disturbance as the air stream velocity increases. However, the 
air shear stress disturbances do not significantly change the wavelength of an ice pattern 
occurring as a result of morphological instability of the ice- water interface. The model herein 
predicts that a centimeter scale ice pattern will appear, and its wavelength will decrease with 
increasing air stream velocity. Moreover, the ice pattern will translate towards the water 
source with time. At higher airspeed, the theoretical predictions obtained here might be 
relevant to the experiments for surface roughness characteristics associated with leading 
edge ice accretion on a NACA 0012 airfoil at a 0-deg angle of attack.— In that experiment, 
the height and spacing of roughness elements were measured with various icing parameters 
in glaze icing conditions. It was observed that the roughness spacing is about 1 mm, and 
that smooth-to-rough zones move upstream towards the stagnation region with time. 

Here, first we mention some differences between previous wet icing models^^"— and the 
current model: (1) The undisturbed part of water film velocity profile derived herein, m;* = 

is the same as that used in the shallow- water icing model.— However, the disturbed part 
of water fiow due to the disturbance of the ice-water interface is taken into account herein. 
(2) The current undisturbed part of temperature in the water film has a linear profile, as 
used in the models.—""— However, in this model, the disturbed part T/ due to the disturbance 
of the ice- water interface is considered, as shown in Eq. ( 150|) . Since the Peclet number herein 
is large, the disturbed part of temperature distribution in the water film is affected by the 
advection due to uu and v'l, as indicated in the terms with Pe/ in Eq. ( l38l) . (3) In previous 
icing models, detailed calculations concerning the effect of the interaction between the air 
and water flows on the temperature distribution in the water fllm were not carried out. The 
air shear stress disturbances influence the disturbed part of stream function of the water 
flow, fi. As a result, the disturbed part of temperature distribution in the water fllm. Hi, is 
affected by both air and water flows. 

If we neglect the disturbed part of temperature distribution in the water fllm flow 
and focus on only the influence of the temperature distribution in the air on the growth 
condition of the ice-water interface disturbance, Eq. ( 1T6|1 may be replaced by L{V + 
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d(/dt) = —KadTa/dy\y=^. Linearizing this equation at y = Jiq, the zeroth order yields 
V = —KaTcx,/{L6o/Ga*), which is identical to Eq. fl32|) . It is found that since the undis- 
turbed part of ice growth rate, V, does not include any parameter associated with the water 
film, V is determined without considering the details of the water film, and the heat trans- 
fer through the air boundary layer is the deciding factor in V. From the first order in ^j., 



dC*/dU = q',, = h'JK yields ai'^ = (KM^^ and Vp, = -ai'^/k^ = -{h'JK)^^ /ku. When 



we neglect the details of the water film, {h'^/h^)^^^ (solid curves) in Fig. |6] (c) directly rep- 
resents the amplification rate. However, there is a significant difference between a^'^ in Fig. 
m (a) and that in Fig. |H] (c). a^^ in Fig. H] (a) takes into account the effect of the disturbed 
part of temperature distribution in the water layer as well as that in the air boundary layer 
on the ice growth conditions. On the other hand, al^'* = {h'^/hx)^'^^ in Fig. IH] (c) is obtained 
without considering the details of the disturbed temperature distribution in the water film. 
This suggests that the wavelength of ice pattern and its translation velocities shown in Ta- 
ble [TTl cannot be evaluated correctly if we neglect the details of the water film. It should be 
emphasized that the same issue arose in a previous paper;— If we neglected the influence 
of the disturbed temperature distribution in the water film flow on the growth condition 
for the ice-water interface disturbance, the amplification rate al had positive values for 
all wave numbers and hence a characteristic wavelength of icicle ripples was not obtained. 
A centimeter scale of ripples in wavelength was obtained from only a^j'^ being taken into 
account the disturbed temperature distribution in the water film.— 

Second, we mention the temperature at the ice-water interface, T,, in Eq. ( lT5l) . Within 
linear stability analysis, Tj = Tsi + AT,;, where Tsi is the temperature at an undisturbed 
ice-water interface and AT^/ is a deviation from it when the ice-water interface is disturbed. 
Its dimensionless form ^T^u = Ini[ATs;/(Ts; — Tj^)] can be expressed as follows by evaluating 
Eq. fl50l) at the disturbed ice- water interface = C*: 



ATsi* = 5h{U) \ {Hf'\y^=o - 1) sm[ku{x* - VpJ^)] + Hl''\y,=QCOs[ku{x^ - VpJ^)] \ . (57) 



Figure [7] (a) represents the isotherms in the water film. The real and imaginary parts of the 
disturbed part of temperature distribution in the water film, and in Eq. ( ISOI) . are 
determined by solving Eq. (l38l) subject to the boundary conditions (jSj) and (H5l) . which 
were derived from the continuity of temperature and heat flux at the water-air interface, 
Eqs. ([T7|l and (fT8|) . respectively. Since the temperature distribution in the water film is 
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affected by botli air and water flows, ATgi* varies. Wlien tlie ice-water interface is flat, 
ATgu must be zero. Indeed, Hj^^'^\y^=Q — )■ 1 and Hl^^\y_^=o — ?■ in tlie limit ki^ — )■ were 
numerically confirmed. 
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FIG. 7. For Q/lu, = 1000 [(ml/h)/cm], Mqo = 20 m/s and = 0.05, (a) represents the isotherms 

in the water film and (b) the isotherms in the ice, for the boundary condition Ti\y=^ = Ta\y=s^ = Tia 
(constant), (c) represents the variation of disturbed heat flux \q[J\, \q'sJ\ and \q'aJ\ against the free 
stream velocity Uoo- For the boundary condition Ti\y=(^ = Ts\y=c_ = Tsi (constant), (d) represents 
the isotherms in the water film and (e) dimensionless amplification rate a^^ versus dimensionless 
wave number ka*. The numbers in the isotherms are the values of dimensionless temperatures in 
water ([50]) and in ice (fST]) . 
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From Eqs. (152|) . ( 153|) and ( 15^ . we define the magnitude of g^^, and as follows: 

\ql\ ^ {idHt^/dy.\y,=of + idHlVdy.\y.=ofV^', (58) 
IgLl = K^k4iHt\,=o - If + (^fi.=o)^}^/^ (59) 

\a - {(G^ //^i.=i - G'^y^%=,r + (Gf + //^)|,.=i)n^/^- (eo) 

Since the water flow is affected by the air flow through the air shear stress disturbances, 
and Ig^^l as well as \q'^J depend on the free stream velocity Uoo- Figure [7| (c) represents the 
variation of \q[J, jg^*! Wa*\ against Uoq. These values were estimated for Q/lw = 1000 
[(ml/h)/cm] and x = 0.1 m, and for ka* at which ai^"^ acquires a maximum value for a given 
Uoo- From Eq. f|T6|) . the disturbed part of the Stephan condition is d(^/dt^ = — g^^. 
Therefore, the net heat flux g^^ — g(,^ determines the amplification rate of the ice-water 
interface disturbance. Since \q'^J\ = \h'^/h^\ from Eq. (13^ . the dashed-dotted curve in Fig. 
[7] (c) is the same as the solid curve in Fig. [6] (b). When both ice- water and water-air 
interfaces are flat and the undisturbed temperature in the water film is the linear profile 
Th = y*i the latent heat released at the ice- water interface, LV , the undisturbed part of 
heat flux at the ice-water interface, KiGi and that at the water-air interface, KaGa, must be 
equal. Hence all of the latent heat is released away from the water-air interface through the 
water film. However, in the case of the disturbed interfaces, the disturbed part of heat flux 
at the ice- water interface, |g;^|, is not necessarily equal to that at the water-air interface, 
|g^^|. Indeed, as shown in Fig. [7](c), |g^^| is much smaller than \q'i^,\. The release of latent 
heat at the water-air interface, Ig^^,!, is limited not only by the temperature gradient at 
the water-air interface but also by the shape of the water-air interface. Hence all of the 
latent heat released at the disturbed ice-water interface cannot be removed at the water-air 
interface and most of it is carried by the flow in the water film. 

The effect of morphological instabilities is to increase the surface area of the phase bound- 
ary and hence to enhance the release of latent heat. The latent heat that is not completely 
removed by air and water flows may lead to a local temperature rise in the supercooled 
water and ice locally. Also, the flow in the water film can carries a supercooled water in 
the interior towards the ice-water interface. Hence, in Fig. [7] (a), not only the isotherms 
in the water film are deformed by the advection terms in Eq. (l38l) but also alternating 
patterns of warming and cooling appear in the neighbourhood of the ice-water interface. 
The characteristic time of the deformation associated with the shear rate is 1/ {uia/ho). On 
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the other hand, the thermal diffusion time associated with a wave number ki is l/{nikf). 
For a disturbed ice- water interface with a 1 cm wavelength shown in Figs. [7] (a) and (d), 
the condition l/iuia/ho) ^ l/{Kikf) is satisfied. Hence, the temperature distribution in 
the neighbourhood of the disturbed ice-water interface is deformed by the advection in the 
water film. In order to avoid a temperature discontinuity, AT^/*, between water and ice, 
the disturbed heat flux q'^^ due to thermal diffusion occurs in the ice. As a result. Fig. [7] 
(b) also shows alternating patterns of warming and cooling in the ice. From the comparison 
of Eqs. ( 157|) and (153|) as well as that of Figs. [7] (b) and [5] (b), if ATsu > (warming), 
then q'g^ < 0, hence the direction of q'^^ is from the ice- water interface to the ice. On the 
other hand, if ATgi^ < (cooling), then q'^^ > 0, hence the direction of q'^^ is from the ice 
to the ice-water interface. However, it should be noted that this disturbed part of heat 
flux in the ice exists only in the vicinity of the ice-water interface, as Eq. fl5T]) shows that 
the disturbed temperature in the ice is exponentially attenuated, and the ice temperature 
approaches = Tgi {Tgi =0 °C for pure water) far from the ice-water interface, as shown in 
Fig. [71(b). 

On the other hand, the isotherms in Fig. [7] (d) are determined from different solutions 
H^''^ and which are obtained by solving Eq. ( l38i) subject to the boundary conditions 
Hi\y^=Q = 1 and Eq. fH5l) . Hi\y^=o = 1 is derived from the condition Ti\y=(^ = Ts|y=^ = Tgi 
(constant). In this case, the temperature at the water-air interface is an unknown to be 
determined, which is deviated from Tia in Eq. f|T7j) and hence the first boundary condition 
in Eq. (15^ is replaced by 

However, since Ka/Ki <^ 1, Ha\r^=Q ~ 1 is a good approximation, which is equivalent to 
the condition Ti\y=^ = Ta\y=^ ~ Tia (constant). As shown in Fig. [7] (d), the isotherms are 
slightly deformed from the undisturbed temperature distribution Ti^, = and the variation 
of temperature in the neighbourhood of the ice-water and water-air interfaces are strongly 
affected by the boundary temperatures, T;*|y^=^^ = Ts*|j^^=^^ = 0.0 (constant) and Ti^\y^=^^ = 
Ta*\y^=^, ~ —1.0 (constant). In particular, the isotherms in the neighbourhood of the ice- 
water interface in Fig. [7] (d) are significantly different from those in Fig. [7] (a). Although 
there exists a shear flow in the water film, the temperature distribution is almost symmetric 
around any protruded part of the ice-water interface. The long arrows at the ice-water 
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interface in Fig. [7] (d) show that the temperature gradient is largest at each protruded 
part of the ice-water interface. This promotes the ice growth at the protruded part than 
at the depressed part, always resulting in an unstable growth of the ice- water interface. 
Consequently, the amplification rate a^J^ obtained from the boundary condition T;|y=^ = 
Ts\y=(^ = Tsi (constant) has positive values for all wave numbers, as shown in Fig. [7] (e). 

From l/{uia/ho) = l/{Kikf), the wave number at which the two time scale equals is given 
by ki^ = y/Pei ~ 20 for Q/lw = 1000 [(ml/h)/cm]. The corresponding wavelength is 150 fim 
from ku = kJiQ = 20 for Uoo = 20 m/s. The effect of the water flow on the isotherms in such 
a microscopic length scale is negligible. Instead, taking into account the Gibbs-Thomson 
effect, the temperature at the ice-water interface is expressed as Tj = Tsi + {TgiT / L)d'^ ( / dx^ , 
from which Hi\y^=Q = 1 — ((io/^o)(^th/^o)^E, is derived. Here Ith = K'i/V and do = TgiTCpi/L^ 
are a macroscopic and microscopic characteristic length, respectively, F is the ice-water 
interface tension and Cpi is the specific heat at constant pressure of the water. Neglecting the 
advection terms in Eq. ( l38l) and solving it subject to the boundary conditions Hi\y^=Q = 1 — 
{do /ho) {Ith/ ho) kf^ and Eq. ( 145|) with fi^\y^^i ^ 0, from Eq. ( 146|) we obtain the amplification 
rate, ai^^ = ki^,{l — {do/ ho){lth/ ho)ki^{l + K[)} , which is just the result of the MuUins-Sekerka 
instability.— Here the condition //*|y,=i ~ means that the water-air interface is nearly fiat 
because the surface tension is so dominant in the microscopic length scale of dendrite that 
the effect of disturbance concerning the dendrite spacing on the shape of the water-air 
interface is negligible, ai^^ acquires a maximum value at k^ = [h^/ {3lthdo{^ + K^)}]^^"^ ~ 10 
for /iq ~ 5 X 10~^ m, /th ~ 10^^ m, do ~ 10~^ m and i^f = 3.92. The dependence of the 
microscopic wavelength predicted from the MuUins-Sekerka instability on the free stream 
velocity Uoo is Amicro = 27r{3/th<^o(l + ^i)V^'^ u^^^ because /th ~ 5o ~ u^^'^ in the model 
herein. This result is contrast to the dependence of the macroscopic wavelength Amacro on 
Moo- Figure |4](b) shows that Amacro oc u'^^^^. It should be noted that the mechanism of the 
macro-scale morphological instability under a supercooled liquid film herein is quite different 
from the dendritic growth. Amacro depends on the water film thickness ho as indicated in 
Tables HI and IIT| while Amicro does not depend on it. 

The use of the boundary condition Ti\y=c = ^s|j/=c = '^si (constant) caused a serious 
problem in a model for the icicle ripple formation.— This condition was used in the model 
proposed by Ogawa and FurukawaM when determining Hi . However, our numerical analysis 



did not reproduce their amplification rate shown in Fig. 4 of Ref. 
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Instead, the numer- 



ically obtained amplification rate had positive values for all wave numbers and there was 



no well-defined maximum amplification rate (for details, see Figure 5 (c) in Ref. llSl ). The 
same issue has already arisen in aircraft icing problems. Tsao and Rothmayer developed a 
physical model to describe the aero-hydro-thermo-dynamic interaction of a cold air bound- 
ary layer with glaze ice sheets and water films.— However, their stability analysis showed 
that the ice-water interface disturbance became unstable for all modes (see Figure 7 in Ref. 



30l ). which indicated that there is no dominant amplification rate to select a preferred wave- 
length. The assumption that the disturbed ice-water interface is at the equilibrium freezing 
temperature was used in their model too. To overcome this issue, the Gibbs-Thomson effect 
was introduced to stabilize the smallest scale icing disturbances.— However, the length scale 
predicted by their theory was much smaller than the ice roughness spacing of the order 
of millimeters observed in the experiments by Shin.— On the other hand, the condition 
Ti\y=(^ = Ts\y=(^ = Tsi (constant) was not used in the model for the icicle ripple formation 
proposed by Ueno^^— when determining Hi. That model was able to predict a centimeter- 
scale wavelength of icicle ripples and upward ripple migration due to an asymmetry in the 
temperature distribution between the upstream and downstream side of any protruded part 
of the ice-water interface, which were confirmed by the experiments.—"^ 

Third, we mention the effect of heat conduction into a substrate beneath an ice sheet 
of finite thickness on the morphological instability. In the model herein, it was assumed 
that the ice region is semi-infinite and the undisturbed part of temperature gradient in the 
ice does not exist. We relax this assumption by including heat conduction into a planer 
aluminum substrate beneath the ice sheet. In Fig. |H] (a), 60 and /sub are the thickness of ice 
and aluminum plate, respectively, T^ub is the temperature between the ice and aluminum 
plate and Tgubo is the temperature of other side of surface of the aluminum plate. Then, the 
undisturbed temperature gradient in the ice is Gs = {Tsi — Tsuh)/bo, and Eq. (146|) is replaced 
by 

(r) dHl"^ COsh{kubo/ho) ( , rr{r)\ A /rqN 

dy^ j/.=o smh(/c;^,Oo//io) ^ ^ 

where Gf = Gg/Gi is the ratio of the undisturbed temperature gradient at the ice-water 

interface in ice to that in water. The ice thickness 60 is determined by integrating the 



following equation subject to an initial condition of 60 = at t = 0> 



.18 



'dt~T bo 17 bo/Ga. ■ ^ ^ 
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FIG. 8. (a) Schematic of an ice growth on an aluminum plate under air and water flows, (b) 
For = 1000 [(ml/h)/cm] and — 20 m/s, dimensionless amplification rate cri ^ versus 

dimensionless wave number /cq* for various ice thickness 6o- The thickness of air boundary layer in 
(a) is defined hy 5 = So/Ga*, where Ga* = 0.413 and 5o = 361 fiin for Uoo — 

20 m/s. 

If other side of surface of the aluminum plate is exposed to ambient cold air, i.e. assuming 
^subo = ^oo, G^i can be expressed asr^ 

Kabof \ Ksub bo J 

where -ft'sub = 237 J/(mKs) is the thermal conductivity of the aluminum plate. 

Figure [S] (b) shows the dimensionless amplification rate ai^^ versus dimensionless wave 
number /c^^, for various ice thickness bo. In the case of b^ = 6o, ai^^ is negative for all wave 
numbers, which means that the ice-water interface disturbances diminish with time. On the 
other hand, in the case of bo = 106o, the ice-water interface disturbances in a finite range 
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of wave numbers become unstable and cr; acquires a maximum value at a wave number. 
Noting that in Eq. (!64|) is zero in the limit Bq ^ and Jiq is the same order as 60 as 
indicated in Table HI when ^ 60 Eq. fl62|) reduces to Eq. fj46l) and the solid curve in 
Fig. M (b) is the same as that for Uoo = 20 m/s in Fig. H] (a). When the ice thickness is 
small during the ice growth, the morphological instability of the ice surface is suppressed 
because the removal of the latent heat due to the conduction into the aluminum plate is 
dominant than that due to the convection by air and water flows. Even in the presence of 
the undisturbed part of temperature gradient Gg in ice, the morphological instability occurs 
when the ice thickness exceeds a critical value. 

Finally, we mention some limitations of the proposed model. First, freshwater icing 
sponginess containing non-negligible amount of liquid water is observed in icicles^iii^ and 
aufeis.— The spongy icing phenomenon is also well-known to the in-fiight icing community.—"^ 
In the experiments of icing wind tunnel using a NACA 0012 airfoil, a considerable variation 
in sponginess (or liquid fraction) with air temperature, wind speed and liquid water content 
was observed.— The model herein cannot explain these experimental results. Therefore, 
it needs to modify an air-water-ice multi-phase system in Fig. [1] to an air-water-spongy 
ice multi-phase system in Fig. [9l where a spongy layer is introduced in between a fully 
water region and a fully ice region. However, since the water region is thin liquid film 
and the latent heat transfer is strongly affected by the existence of the water-air interface, 
the configuration in Fig. [9] is fundamentally different from the mathematical models of 
flow-induced morphological instability of mushy layers developed in previous studies,—"— 
where the liquid region was taken to be semi-infinite. 

The conventional Stefan problem cannot describe the pattern formation observed in na- 
ture, because the dimensional information needed to set the scale of a crystal growth is 
absent.— In other words, if the temperature at the ice-water interface is kept at °C and 
neglecting surface energy, the morphological instabilities occur on arbitrarily small length 
scales given any amount of supercooling.— In practice, surface energy limits the instability 
at some scale and an ice surface under a supercooled water film results in dendritic growth. 
As a result, there can be a possibility of spongy ice formation, in which a portion of the 
surface liquid is incorporated into the dendritic ice matrix.— Hence the spongy layer is a 
mixture of ice and water and its temperature is a thermodynamic equilibrium one.- The 
use of this local equilibrium assumption is appropriate in the interior of the spongy layer 
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at some distance from the tips of dendrites, where the increase in specific surface area of 

micro-scale phase boundaries promotes the release of latent heat into the interstices (pores) 

and hence the level of non-equilibrium can be kept very small.— 

) 
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FIG. 9. Schematic of air-water-spongy ice multi-phase system. 

In Fig. m since the interface between the spongy ice region and water region does not 
have a well defined position on the micro scale of dendritic growth, the spongy ice-water 
interface is defined as the envelope (suitably smoothed) of the dendrite ice matrix.— If 
the ice-water interface in Figs. [7| (a) and (b) is replaced by the spongy ice-water interface, 
Tj in Eq. ( fT5l) is interpreted as the temperature at the spongy ice-water interface. The 
local equilibrium assumption is likely to break down in the neighbourhood of the spongy 
ice-water interface, because if the spongy layer is postulated as a permeable material,— a 
significant effect of the tangential and normal air shear stress disturbances found herein on 
the spongy layer through the thin water film can be expected and the latent heat advected 
by the water flow penetrates into near the spongy ice- water interface. The level of this 
non-equilibrium effect in the macro-scale is negligible at some distance from the spongy 
ice- water interface, where the temperature is °C, as shown in Fig. [7] (b). In order to 
gain a clear understanding of these viewpoints, it is necessary to add equations governing 
local liquid fraction and the internal evolution of the spongy layer to the model herein and 
the effects of the shear stress at a disturbed spongy ice-water interface on the distribution 
of liquid fraction, permeability and penetration depth should be investigated along with a 
study on non-equilibrium coexistence of crystal and shearing liquid flow.— Furthermore, the 
dependence of the liquid fraction on the icing parameter such as air temperature, wind speed 
and liquid water content in the experiments by Lozowski et al.— must be explained. 
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Second, in the linear stability analysis, the amplitude of the ice-water interface distur- 
bance of the most unstable mode increases with time: = exp{aima.xt*)Sb, which affects 
the magnitude of h'^/hx in Eq. (IHH]) . In order to evaluate the value of crimax, it is necessary 
to determine the disturbed flow and temperature fields in the water film which are influ- 
enced by surrounding airflow and temperature fields. However, the linear theory is unable 
to clarify further features related to the development of disturbance. We have to generalize 
the equation d5b{t^)/dt^ = crlmax5fe(t*) to a nonlinear amplitude evolution equation. 

Third, the magnitude of h'^/h^ depends on the shape of the ice-water interface. In 
the normal mode analysis presented here, the values of h'^/hx depend on the supposed 
sinusoidal form in Eq. fl55|) . It is necessary to extend the current model to the problem of 
ice morphology produced on an arbitrary three-dimensional surface such as aircraft wings 
and overhead line cables,— taking into account water flow driven by both gravity and wind 
drag simultaneously. Removing restrictions mentioned above and further extension of the 
current model to practical icing problems will be discussed in later papers. 
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